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One problem which plagues the numerical evaluation of one-loop Feynman diagrams using recursive integration 
by part relations is a numerical instability near exceptional momentum configurations. In this contribution we 
will discuss a generic solution to this problem. As an example we consider the case of forward light-by-light 
scattering. 



1. Introduction 

By using recursion relations based on the inte- 
gration by part techniques [T] a general algorithm 
can be constructed to numerically evaluate the fi- 
nite part of one loop tensor iV-point integrals . 
The numerical implementation is straightforward. 
However, for specific configurations of the exter- 
nal momenta, the recursion relations build up a 
numerical instability. Phcnomcnologically, these 
exceptional kinematics are characterized by such 
configurations as mass thresholds, forward scat- 
tering and planar event where a parameter tends 
to zero (for example the scattering angle) . In this 
contribution, we will discuss how to modify the 
recursion relations to produce expansion formu- 
lae in the small parameter. This is an extension 
of the techniques developed in ref. [H] (see also 
ref. and allows the evaluation of tensor loop 
integrals at/near exceptional momenta configu- 
rations. As an example, we work out in some 
detail forward light by light scattering through a 
massless fermion loop: 77 — > 77. The general- 
ization to massive particles with mass thresholds 
and processes with larger multiplicities is in prin- 
ciple straightforward. 

2. Exceptional Recursion Relations 

The dimensionality of space-time fixes the 
number of momentum vectors needed to form a 
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basis set in Minkowski space. As a consequence 
the integration by part recursion relations fall into 
three categories depending on the number of ex- 
ternal legs, N, in the loop graph. These cate- 
gories are: N < 5, N = 6 and N > 7. 

At the center of the recursion relations is the 
so-called kinematic matrix 

Sij = (ft - qj) 2 - rn\ - m] , (1) 

where qi is the momentum flow through and rrii 
the mass of propagator i. For the three afore- 
mentioned categories of recursion relations, the 
kinematic matrix has the following properties: 

N < 5: Both the inverse of the kinematic matrix 
and the Gram determinant exist for non- 
exceptional momentum configurations. The 
Gram determinant is proportional to 

5 = E^ = E^- ( 2 ) 

i ij 

For momentum configurations close to ex- 
ceptional kinematics the B parameter tends 
to zero. Both the recursion relations and 
the basis set of loop integrals become nu- 
merically unstable in this limit. We need to 
setup a recursion relation which generates 
an expansion in B. 

N = 6: In this case, the inverse of the kinematic 
matrix still exists for non-exceptional mo- 
mentum configurations. However, B = 
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^ S^ 1 = for all momentum configura- 
tions. This leads to a different set of re- 
cursion relations. Momentum configura- 
tions close to exceptional kinematics are 
now characterized by small eigenvalues of 
the kinematic matrix. In other words, the 
numerical determination of its inverse be- 
comes unstable. We need to setup a recur- 
sion relation which generates an expansion 
in the small eigenvalue(s) of the kinematic 
matrix. 

N > 7: Here the kinematic matrix is singular for 
all momentum configurations. Note that we 
do not need any special recursion relations 
since the standard recursion relations are 
numerically stable. 

As we look at changing kinematic behaviour 
as the multiplicity increases, we also see how 
to construct the recursion relations near excep- 
tional kinematics. Specifically, when we are ex- 
actly at an exceptional momentum configuration 
for N < 5, i.e. B = 0, we can use the N = 6 re- 
cursion relations which were constructed on the 
premise that B = 0. Similarly, for N = 6, at the 
exceptional momentum configurations the N > 7 
recursion relations can be used. 

However, we are interested in the phase space 
regions close to the exceptional momentum con- 
figurations. It is clear how to construct the ap- 
propriate recursion relations in these regions. We 
need to rewrite the N < 5 recursion relation as 
the N = 6 recursion relation plus terms propor- 
tional to B (i.e. the Gram determinant). As we 
will see in the next section, this leads to an ex- 
pansion formula in B, enabling us to evaluate 
the loop graph near the exceptional configura- 
tions with arbitrary precision. Equivalently, we 
can rewrite the N = 6 recursion relation into the 
N > 7 recursion relations plus terms proportional 
to the small eigenvalue (s). This leads to an ex- 
pansion in the small eigenvalue(s). 

3. Forward 77 — > 77 Scattering 

In this section we work out the simple exam- 
ple of light-by-light scattering through a mass- 
less fermion loop. The analytic answers for the 



helicity amplitudes are well known. Some helic- 
ity amplitudes, such as 7+7+ — > 7^7^ are con- 
stants, independent of the underlying kinemat- 
ics. Other helicity configurations give simple be- 
haviour which is logarithmically divergent as the 
scattering angle goes to zero. Applying the recur- 
sion algorithm of Ref. [2] to this process leads to a 
numerical evaluation of the amplitudes. However, 
when the scattering angle becomes small the algo- 
rithm becomes numerically unstable (see Figs. 4 
and 5). The reason is twofold. First of all some 
of the recursion relations depend on the inverse 
of the B parameter 
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which for small scattering angle 8 tends to zero. 
Secondly, one of the basis integrals, i.e. the end- 
point of the recursion relation, becomes unstable. 
More precisely, the six-dimensional box integral 
has a numerically unstable form (see Figs. 2 and 
3): 

lim 7(6; 1,1, 1,1) 

sB— >0 

= lim 

sB^O tsB 

in 
s 

As discussed in Sec. 2, the solution is to derive an 
exceptional recursion relation: the B-expansion 
relations. This is quite straightforward. We use 
the N = 6 recursion relation plus a correction 
term proportional to B: 

I{D;v 1 ,v 2 ^z,Vi) = 

biI(D;vi - 1,^2,^3,^4) 
+ b 2 I{D;v 1 ,is 2 - 1,^3, ^4) 

+ &4/(Z?;z/i,Z/ 2 ,f3, ^4 - 1) 

+ B(D + l-a)I{D + 2;is 1 ,v 2 ,V3,V4) (5) 
where 
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ate propagator. This integral depends on the two 
Mandelstam invariants, s and t. 
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This relation is exact for N < 6 in any momen- 
tum configuration. However, unlike the usual re- 
cursion relations it has no termination point. For 
example 

7(6; 1,1, 1,1) = 

6i/(6 - 2e; 0, 1, 1, 1) + b 2 I{<6 - 2e; 1,0, 1, 1) 
+ 6 3 7(6 - 2e; 1, 1, 0, 1) + b 4 I(6 - 2e; 1, 1, 1, 0) 
+ (3-2e)B7(8-2e;l,l,l,l) . (7) 

The box integral 7(8 — 2e; 1,1, 1,1) can be ex- 
pressed in terms of 8-dimcnsional triangles (pro- 
portional to B) and a 10-dimcnsional box integral 
(proportional to B 2 ). By iterating this process, 
the scalar box is expressed as a series of triangles 
with increasing factors of B. In the case of the 
massless 6-dimensional box we simply get: 
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Y B m a m I(6 + 2m - 2e; 1, 0, 1, 1) 

_m=0 

+ 0(B M+1 ) 



(8) 



with a m the expansion coefficients generated from 
the last term in eq. 7. In other words we can eval- 
uate the 6-dimensional box integral in the forward 
scattering region as an series expansion in B, with 
each coefficient given as triangle integrals with 
one offshell leg. If the scattering angle is small, 
not many iterations are needed for an accurate 
evaluation of the box diagram. This is illustrated 
in Fig. 1. Of course, for large B, the unmodi- 
fied integration by parts relations are completely 
stable. 

Note that the triangles appearing in the B- 
expansion relations are all UV divergent. How- 
ever, order by order in B these divergences cancel 
in the sums (because by construction the original 
integrals are UV finite.) This means wc can use 
a regulated expressions for a m and the triangles 
integrals. A complication is that the factor a m 
itself is dependent on the dimension. This means 
that the limit e — > has to be taken after the UV 
divergences have been cancelled by adding differ- 
ent triangles. An important property of the coef- 
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Figure 1. The value of — sB vs the depth of ex- 
pansion m needed to evaluate the scalar box in- 
tegral 7(6; 1,1, 1,1) with a relative accuracy of 
10" 25 , 1(T 15 and 1CT 5 . 



ficient a m is that it is a function of the expansion 
depth m and the value of D — a where both D 
and o are given by the triangle the coefficient is 
multiplying. The regulated expression to be used 
in the 7?-expansion relations is given by 

J m {D;vi,v 2 ,vo > ) = a m I{D + 2m\vi,V2,vz) = 

2 m (-l) g Q 2 J ) / ) U/2 + n-m) 
r(cr + 2n) V JnK JnK ' > m 

x (2*(cr + 2n) + *(n + 1) + *(cr/2 + n-m) 

-*(i/ 2 + n) - *(za. + n) - (ct/2 + n) 

-\og(-Q 2 )+ lE ) (9) 

where n is given by the relation D = 2(a + n), the 
Pochhammer's symbol {n) = T(n + m)/Y(n), 
the Eulcr constant and the Digamma function 
ty(x) = d log (r(x)) jdx. As indicated in the no- 
tation, in the algorithmic expansion we only need 
to keep track of the tuple (m, D, v\, 1^2,^3) in or- 
der to evaluate J m (D; v\, v 2l ^3). For example, 
the B-expansion of 7(8; 2, 1,1,1) using Eq. (5) 
now becomes 

M 

7(8, 2, 1, 1, 1) = 6!^S m J m (8; 0, 1, 1, 1) 



m=0 



M 



b 2 J2B m [J m (8; 2, 0, 1, l)+6i J m (8; 1, 0, 1, 1) 



m=0 



4 




1e-12 




1e-13 



Standard recursion 
B expansion 



1e-9 



1e-5 



1e-1 



-B s 



Figure 2. The value of —sB vs the ratio of the 
standard recursion result and the S-expansion re- 
sult for 7(6; 1,1, 1,1). 



Figure 3. The value of —sB vs the imaginary part 
of the 6-dimensional box integral, 1(6, 1, 1, 1, 1). 
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+ b 3 J2B m [J m (S; 2, 1, 0, l)+6iJ m (8; 1, 1, 0, 1)] 



m=0 
M 



+ b 4 J2B m [J m (8; 2, 1, 1, 0)+6i J m (8; 1, 1, 1, 0)] 

m=0 

+ 0(B M+1 ) . (10) 



Note that for N = 5 the expansion becomes a bit 
more involved, but is still rather straightforward. 

By looking at the ratio of the standard re- 
cursion relation and the £?-expansion relation in 
Fig. 2 we can see clearly the numerical insta- 
bility of the standard recursion relation when 



-sB < 10" 



As can be seen from Fig. 1 we 



already can neglect terms of order B 3 in the 
B-cxpansion relation to achieve an accuracy of 
10 -25 in the integral evaluation. Also shown, in 
Fig. 3 is the imaginary part of the scalar integral 
using both recursion relations. As can be seen, 
the B-expansion relation can be used over a large 
range of values provided the expansion depth is 
sufficient (as indicited by Fig. 1). 

We now can evaluate the 77 — > 77 cross section 
in the forward region using both the standard and 
5-cxpansion relations. We define the normalized 



amplitude A(6) 



M(6) = —A(9) 



(11) 



where e is the coupling strength of the photon- 
fcrmion vertex. 

In Fig. 4 we show the helicity amplitude for 
7~7~ scattering. In this case we have 



A++--— (8) = 1 . 



(12) 



As can be seen in Fig. 4, the standard recursion 
relation has a sudden loss in accuracy at around 
6 w 3 x 10~ 3 . On the other hand, the B ex- 
pansion formula gives the correct result to within 
arbitrary precision. 

Similarly, in Fig. 5 we show the helicity scat- 
tering 7+7+ — > 7+7+. The normalized helicity 
amplitude is now more complicated 
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where c = cos 6. The quantitative behaviour is 
identical to the previous case; the B improved 
recursion relations give an accurate result in the 
domain where numerical instabilities render the 
standard approach invalid. 
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Figure 4. The scattering angle vs the normalized 
amplitude for 7+7+ — * 7^7^. 



Figure 5. The scattering angle vs the normalized 
amplitude for 7+7+ — > 7+7+. 



In summary, the i?-expansion method can be 
used to extrapolate the algorithmic integration by 
parts method to all values of 9 without any loss 
of accuracy. 

4. Conclusions 

The method outlined can readily be generalized 
to more complicated processes. As explained in 
sec. 2 this is done by rewriting the recursion re- 
lations as expansion in the small B-parameter or 
small eigenvalues of kinematic matrix. Because 
the recursion relations are numerically solved for 
each phase space point, the method allows ex- 
ceptional kinematic configuration to be automat- 
ically detected by the algorithm on an event-by- 
event basis. If needed the algorithm can decide 
to use the expansion method without any knowl- 
edge of the underlying physics (e.g. threshold 
region, planar event configuration or other more 
complicated configurations). This is highly de- 
sirable because of the complexity of processes 
we are ultimately interested in. For example, 
four quark final states at hadron colliders such 
as PP — > ti + bb have multiple mass thresholds 
and numerous other kinematic exceptional con- 
figurations. Similarly processes like PP — > W + 4 
jets requiring the evaluation of the complicated 
one-loop diagrams involving six partons plus a 



vector boson for which the exceptional kinematic 
configurations are difficult to comprehend. 

The algorithm outlined in these proceedings 
augments the integration by parts algorithm of 
Ref. 2 j. The combined numerical procedure 
should be able to calculate arbitrarily compli- 
cated one-loop amplitudes in all regions of phase 
space with arbitrary precision. 
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